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1. ABSTRACT 


A self-adaptive grid procedure for efficient computarion of three-dimensional complex flow fields is 
described. The method is based on variational principles to minimize the energy of a spring system anal- 
ogy which redistributes the grid points. Grid control parameters are determined by specifying maximum 
and minimum grid spacings. Multidirectional adaptation is achieved by splitting the procedure into a 
sequence of successive applications of a unidirectional adaptation. One-sided, two-directional constraints 
for orthogonality and smoothness are used to enhance the efficiency of the method. Feasibility of the 
scheme is demonstrated by application to a multinozzle, afterbody, plume flow field. Application of the 
algorithm for initial grid generation is illustrated by constructing a three-dimensional grid about a "bump- 
like" geometry. 


2. INTRODUCTION 


The effective design of future aircraft, rockets, and aerospace vehicles which fly at high speeds 
requires realistic consideration of three-dimensional (3-D) complex flow fields. Analysis of these prob- 
lems by means of simulation techniques demands enhanced grid distribution methodologies to obtain effi- 
cient and accurate resolution of rapidly varying solutions to the pertinent partial differential equations 
(PDE). 

Studies of solution-adaptive grid methods are significant because of their potential for improving the 
efficiency and accuracy of the numerical methods and for reducing computer memory requirements. Con- 
tinuing advances in solution algorithms and computer power have considerably improved the efficiency of 
the solution, but the accuracy of the solution still depends directly on the distribution of grid points over 
the physical space in question. 

In recent years, a considerable amount of literature has been devoted to the subject of adaptive 
methods. General surveys by Anderson (ref. 1) in 1983, by Babuska et al. (ref. 2) in 1983, by 
Thompson (ref. 3) in 1985, and by Eiseman (ref. 4) in 1987 summarizes much of this activity. One of the 
important tools which has been used to cluster the grid points in the physical space and which has been 
reviewed in detail, in references 1 and 3, is to redistribute the grid points along a fixed line so that some 
positive weight function is equidistributed over the line. Similar adaptation along a family of fixed lines 
can be applied to multidimensional adaptation (ref. 5). Without additional constraints these methods can 
produce undesirable results in terms of grid skewness. By introducing an orthogonality constraint, 
Anderson (ref. 6) and Dwyer (ref. 7) have been able to control this problem to some extent in two- 
dimensional (2-D) applications; Eiseman (ref. 8) and Nakahashi and Deiwert (refs. 9,10) have introduced 


’National Research Council Research Associate, Member AIAA. 
Thief, Aerothermodynamic Branch, Associate Fellow. 


1 


a practical multidimensional adaptation using successive applications of a one-dimensional (1-D) adapta- 
tion separately in each of the curvilinear coordinate directions. In the latter a torsion-tension spring system 
analogy is used incorporating resistance to movement away from orthogonality. The orthogonality and 
smoothness of the above spring system analogy have been controlled by several incorporated parameters 
and constants whose values are determined by the effort of the user through numerical experiment. To 
diminish the degree of empiricism intrinsic to this approach, particularly in 3-D, Nakahashi and Deiwert 
(ref. 1 1) have suggested some new measures to relate the value of the parameters directly to the solution 
data. This new improvement has significantly increased the accuracy and efficiency in the study of 2-D 
airfoil problems. 

This paper presents a practical 3-D solution-adaptive algorithm for the redistribution of grid points to 
optimal positions and is suitable for modem computer architectures. This work extends the 2-D self- 
adaptive approach of reference 1 1 to 3-D. The control of grid skewness is achieved by one-sided torsion 
and control of spacing by two-sided tension springs. The procedure can be used independently of the 
PDE-algorithm and/or can be coupled with a variety of existing PDE-solvers, both for steady and unsteady 
cases. It is suitable for use with zonal approach. The multidirectional adaptation is achieved using the 
concept of splitting (or fractional steps) in all aspects of the approach. 

A user friendly 3-D code which implements this algorithm has been developed. The code automati- 
cally calculates most of the incorporated parameters, which affects orthogonality and smoothness via user- 
specified constants denoting the maximum and minimum grid spacing. The user selects direction of 
adaptation along a curvilinear grid line in the physical space. Grid points are then redistributed to their 
optimal positions along a family of these lines whose collection sweeps a user- specified curvilinear grid 
plane. The procedure covers the entire physical space along each of these planes in a given marching 
direction. The feasibility of the method as a solution-adaptive scheme is illustrated by its application to a 
realistic three-dimensional problem concerning a generic two-nozzle afterbody plume flow. The outer 
flow field is complicated by the three-dimensional body shape interacting with multiple jet streams, three- 
dimensional barrel shocks, Mach disc, mixing shear layers and recompression shocks. Application of the 
scheme as an initial grid generation is illustrated by generating a 3-D grid about a bump. 


3. METHODOLOGY 


Here we briefly review the conceptual approach, based on variational principles that are the main 
ingredient of the adaptive grid scheme. Then a more relaxed version, although not in content, of some of 
the mathematical constraints are considered to produce an efficient and practical solution-adaptive 
algorithm. r 

The problem of grid genertion and adaptation is to establish the mapping relation x(^) or £(x) 
between the physical and computational space. An optimum distribution is one in which the solution error 
or some approximate measure of this error, here denoted by a positive weight function w(x), is uniformly 
distributed over all grid points. The solution data should be adjusted on the new grid either by interpola- 
non or incorporation of the grid speeds in conjunction with the transformation of the time derivatives. 

This concept is discussed in reference 1 in detail; here we will use the following relations taken from this 
reference for our immediate purpose: 
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discrete case; 


( 1 ) 


r*** 

w(x)dx = constant, 

» X 


AxjWj = constant 


2 

x^w(^) = constant, Ij = w(£)x^ d£ ( 2 ) 

The initial part of equation (2) is the Euler-Lagrange equation for minimization of the second integral there. 
It can also be interpreted as minimization of the energy of a system of springs with constant w(c) between 
each pair. We also have 




w {%)) 


(3) 


where N is the total number of grid intervals. 

The weighting function w(x) is most commonly experssed as a linear combination of positive func- 
tions, Mi's, related to the first or higher gradient of the solution, 

w(x) = M 0 + ciMi + c 2 M 2 + ... + c m M m , Ci > 0 (4) 

See also suggestions by White (ref. 12). 

In multiple dimensions, the adaptation in different directions should be coupled to maintain sufficient 
smoothness in the grid. Brackbill and Saltzman (refs. 13-15) in their variational approach minimize a lin- 
ear combination of integrals where each is a measure of some grid properties; their formulation is given 
by: 


Io = 



1 = 

= x v 


4- ?l s I s + ^-oIq 

(5) 

I w = Jw(x)J dx 

m 

[ w(^)J 2 


( 6 ) 

( 3 1 



f 3 ^ 




dx = 

■J 

X J ii 

Jd^ 

(7) 

J V i=l J 


Vi=l ) 



+ (J 23 ) + (J 31 ) 

]dx 

— 1 

• 

[[(J 12 ) 2 + (J 23 ) 2 + (J 31 ) 2 ]J 2 d^ 

( 8 ) 

Jy = v? ■ 

* 

j = 

(det | Jjj 

il ) 1/2 

(9) 


J is the Jacobian of the transformation, A. w , and Xq are specified coefficients, dx and d^ are differen- 
tial volumes in physical and computational space, respectively. Larger values of these coefficients put 
more emphasis upon their corresponding weighted volume integrals. The I w forces the cell size to be 
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small where w(x) is large and vice versa; I s reinforces the smoothness of the grid; and finally, I Q con- 
strains the orthogonal behavior of the grid. 

Minimization of I yields the Euler-Lagrange equations and can be obtained by application of the 
following operators to the integrand of I: 


E‘ 




i= 1,2,3 


( 10 ) 


Equation (10), which establishes the mapping relation, would now be a system of three coupled elliptic 
PDEs whose solutions would be responsible for generating smooth, clustered or stretched, and orthogonal 
curvilinear coordinates in the physical space. The solution of the previous equations in practice is in itself 
a rather complicated task. Here we will show how a spring system analogy whose constants are functions 
of the solution itself can be used to relax multidimensional elliptic equation (10). The procedure illustrated 
is an extension of the self-adaptive method of Nakahashi and Deiwert (ref. 9) to 3-D problems. 


4. UNIDIRECTIONAL SELF-ADAPTIVE SCHEME 


In this section we will derive the basic equations for grid adaptation. For simplicity of illustration, 
consider a 3-D flow field with adaptation along the curvilinear coordinate £. Also for clarity of our nota- 
tion assume that indices i, j, and k correspond to the rj, and C, coordinates in the computational space 
and let ^ denote the curvilinear plane & = constant, which is a transformation of the \ = i of the com- 
putational plane. The denotes the curvilinear line which is the intersection of the ru and ^-planes. 
Finally, let the point p(i,j,k) of computational space correspond to the point (x; ; k, Yi i k z; ; t ) of the 
physical space and let fy >k be the value of a function f at p. J d * 

In figure 1, grid points on the line are free to move along a line whose configuration is fixed. 
Each point A on ^k is suspended by two tension springs which connect A to points B and C and 
whose spring constants (weighting functions) are wj-ij* and wy* and w is 


w i,j,k = M 0 + rMf , 


Mj = 


fj,j,k ~ *MIN 
*MAX - ffvHN 


(ID 


where fMiN and fMAX are the minimum and maximum of nonnegative driving function f with its dis- 
crete values denoted by fy )k ; r and b are positive constants; and Mo is used to maintain grid smoothness 
and is equal to unity in most cases. The function f > 0 is related to the gradient of pressure, density, 
temperature, Mach number or other vital solution variables and/or a linear combination of these quantities. 
This choice of f is normally used to define the stretching and clustering of grid points, for other pur- 
poses, such as distributions across wall-bounded shear layers a geometric (exponential) function and 
around curved surfaces, the local curvature would be used. 

constants r and b in equation (1 1) are critical in obtaining the degree of grid clustering along the 
grid lines. Determining these constants by means of experiment, particularly in multidimensional 
application, is a tedious task. Here an automatic method of calculating them is considered. This would be 
one area of computation which makes the scheme "self-adaptive." In this approach (see also ref 1 1) a 
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desired maximum and minimum grid spacing Asmax and Asmin. respectively, are specified by the user; 
then the grid spacing is determined by equation (3) and written for three dimensions as: 


( N 


V 1 


A Si,j, k = L| 


w 


i.j.k 


X (w n,j,k) 1 
n=l J 


( 12 ) 


where sjj k is the approximated arc secant length and calculated from point (i,j,k) along the line £,j >k . 
From equations (1 1) and (12), Max Mj the maximum of Mj is equal to unity and Max wjjjc = 1 + r for 
some i = i*, 1 < i* < N when Mo s 1. The ratio of Max Wjj* to Min Wjj ik; , from equation (12), is 
proportional to the ratio of Max As; j ik to Min Asj jjc- These observations all suggest that r should be 
defined as 


^ s max 

As M1N 


(13) 


For practical convenience Asmax and Asmin are given as multiples of the uniform spacing value L/N. 


The constant b is iteratively found such that equation (12) is satisfied; this means we find some 
b = 5 so that y(5) = Min Asy,* has the value \|/(5) = Asmin- This can be done by application of the 
Newton-Raphson iterative scheme; for the (v + l)th iteration of b we have: 

b<*-l) = b (,) + Ab, Ab=(As MIN -<|/(b))/(|j£) (14) 


ch|/ _ 
3b ~ 


E>Vi + D\|/ 2 


(15) 


Min As: : t u 

1X1 ' r Max Aw::;: at i = i* 


(16) 


N 


( Max W; j k \ v M n logM n . 

0^2 = r I r J (Min As i,j,k> Zj — . at 1 = 1" 

( w n.j,k) 


n=l 


(17) 


In equations (16) and (17) we have used the fact that Min Asijjc = Max wjj ik , at i = i* from equa- 
tion (12). Also from Max Mj = 1 and Max wjj jk = 1 + r we have Dxj/i = 0, and D \|/2 is given by 


N 

d V 2 = r(1 ^ ^ (Min As i j k ) 

n=l 


NtflogM, 

(w i,j.k) 2 


i = i* 


(18) 


A force to control inclinations of the r\ and C, coordinates is given by the torsion springs attached at 
points D and E in figure 1; these torsions enforce the inclinations of lines DA and EA with respect to 
prescribed reference lines DAi and EA 2 . If the coefficients for torsion springs are denoted by Ci and C 2 
a mathematical statement of the force is given by 

^torsion = _ C 1 0 i j_i k - C 2 0y ik _i (19) 
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^ EA rT cdve,y - me «*- 

EAH, respectively. The DA? itself k k • u a ^ ianc ^^2 on the planes DAB and 

maintain smoothness along the r|-line and 2) a line nomalto th^hmof DCfi* S' 'T&f. + 1} on 

EHE, J E= + (i, j n + 1, *e VIi\T Hne ed oIh rage f° f ^ eX , tension of QE and normal to the plane 

could be used as well. These springs affect the nrrhnT !j neS ’ SUch as streamlines or shocks, etc., 

equations (7) and (8) The effects are incomnrat H ogonallty and sm oothness constraints described by 
ellipticity in the t, m 3 ° ne - sided mann ^ without introduction of 

balanc^^ Wi,h neW ^ ° f -J* - — * a 

W >.J,k^ s i+l,j,k - S iJik ) - Wi_ x j k ( Si j k - Si_ 1Jik ) - Cj Gj j_j k - C 2 0i j k _! = 0 (20) 

concept is possible, if they are specified along each Inline as follow* aelf-adaptivtty 


N 


Cl ’ C 2^2 - tt = X W ".J. 


n=l 


( 21 ) 


r* h toe X me d adapSd ^ KsfslsSTo afand It* ™?r l “ f ' ension s P ri "S “"a, ants along the 
fomtly over the flow fidd h) ,ha " (Cl and ° 2) where ,he '««* is fixed uni- 

The inclination angles are approximated by 


0 i,M,k = — 1 L . Sj °' k . 0. ^ sy.k 


DA, * 0 i,j,k-l = 


DA, 


( 22 ) 


respectively. Equation ^O) reducesTo '° ^ intersectlon of reference lines DAj and DA 2 , 
W, ’ J ' kSi+1 ’ J ' k " (Wi+1, i’ k + W -Uk + V,, k + T, j k _,) Si j k + Wi _, j kSi _ l j k 

= -T i.j-l.k s ij,k - T i,j,k-l s i,j, k (23) 

T i.j-I,k = (aX,,)/DA, , T iijik _j = (cd 2 )/DA 2 (24) 

In\^Ts' approach, ^niy^h^t^sfon^wces 6 ^ 3 ^^)!^^ 8 ^’ 11 ^ "* **** - « "e readily solved. 

-ion a, the coordinate. This aliows a 
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eighter rjj or ^k-plane) stepping along the j or the k index. However, if the influence of torsion forces 
from both sides (Ty.ik and Ty+i^k) and (Tyjc-i, Ty t k+i) are considered simultaneously, equation (23) 
is no longer a tridiagonal system of equations. This increases the computational effort without additional 
benefit. 

It should be noted that equation (23) without inclusion of torsion forces is a 1-D elliptic equation, 
similar to equation (2), whose solution Sj j ; k is a monotonically increasing function along the ^j^-line. 
With torsion forces included, however, the equidistribution concept is disturbed. This may disturb the 
monotonicity of the solution when inappropriate values of the constants, and X 2 , and the inclination 
angles are used. As a result the grid points may become extremely clustered, coarse and/or even crossed 
and overlapped. To avoid this, one should monitor the values of syjc and Sj for possible 
overlappings; e.g., Si^k < sj+ij^ for some i = / and Sj < si+ijjc for some i = m. A good policy 
is to modify the values of s'y ,k and s'i jjc by equating them to the average of their local maximum and 
minimum in the overlapping region. Otherwise, in the case of severe clustering or coarseness weighting 
functions wy^ should be relaxed (decreased or increased) in proportion to a measure of the fractions 
Asyjc/AsMlN or As; j^/Asmax* respectively. These considerations are automatically done in any course 
of adaptation; however, in the case where monotonicity of the solution values, sy^, fails, the procedure is 
terminated indicating some smaller values of the constants Xi and X 2 be specified. Equation (23) is 
solved for each sy^, iterating until converged values of syjj to some specified tolerance are obtained. 


4.1 Outline of the Procedure 

This section summarizes the basic steps for grid adaptation along principle adaptation directions; i.e., 
a family ^-coordinate lines lying on the £k-plane stepping along the j -index from j = j s to j=j E and 
then marching along the £k-plane from k = ks to k = kg (definitions of these indices will be given next). 
The family of these adapted planes when marched along the k-index covers the entire field. The summary 
is as follows: 1) specify the rectangular box R of the computational space to be adapted, i.e., the bound- 
ary indices is, iE, js> jE> ks, and ke, where S and E denote starting and ending indices of R; 2) specify 
maximum and minimum allowable grid spacing Asmax Asmin, to control the density of the redistributed 
points; 3) specify the torsion spring coefficients Xj and X% 4) march in the k direction (i.e., on each 
£k-plane) beginning at ks; 5) march in the j-direction (i.e., on each £jk-line lying on the £-plane) 
beginning at js; 6) for each k and j calculate arc secant length sy,k and wy^ (from eq. (14)) for 
is ~ i ~ ie; 7) calculate sy,k and Sjj,k> Ty.i^, Tyk-i and then correct for any overlapping of Si,j,k 
and sj jjc; 8) solve the tridiagonal system of equation (23) for new values of sy^ and monitor their 
monotonicity; if it fails the program is terminated, then return to step (3) and select a new set of parameters 
(A-i and X 2 ) and repeat; 9) interpolate fyjc on syjc from the flow solution data on the old grid points; 
and 10) check the convergence of sy^; if not converged keep returning to step 6. 


4.2 Multidirectional Adaptations 

Multidirectional adaptation is achieved by sequential application procedure into a sequence of unidi- 
rectional adaptations. For instance, the grid points are first adapted to the solution along each 
^-coordinate, followed by adaptation along each Ti-coordinate, and then for each ^-coordinate. This is 
analogous to the fractional step schemes for PDEs. It should be recalled that the earlier concept is main- 
tained because of the consideration of torsion forces from one side only. 
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The selected order of adaptation is arbitrary, but the resulting adapted grid is not necessarily uniaue 
One reason is that equation (23) is nonsymmetric and the choice of the boundary of the region ^andthe' 

torsion parameters and Xi play a role in deciding the final grid. 

w?°°l P £ IlCy f ° r Cho ° sing t , he first adaptation direction would be choosing the coordinate for which 
fL , , Jith flowproperties are large. The solution field is interpolated into a newly adapted grid using 

pffi d t°ri? er ’ 0ne ~ dim ^ nS ! 0nal Lagrange interpolation after each unidirectional adaptation and ismore 
necessa^ f ^^grfd^n tl ^ lmens lona ! Relation. For unsteady problems, no inflation would be 

SeS of .hfsolu&S afKr eaCh adap,a “° n a " d tacMporaKd ■“> 

The 3-D code which implements this algorithm enjoys a modular structure. The user can specify anv 
direction of coordinate line for the adaptation and either of two planes of adaptation as discussed earlier- ' 
then the god and solution data are swapped so that the standard procedure oU he previous section can be’ 
applied This portion of the algorithm is another important user friendly feature of this code which auto 
manually maps the inmal input grid and solution data with user-specified directions of aS^rion into 
y or ered data. The adaptation can then proceed along the principal adaptation directions and that 
after completion the data are returned to their original order before output, by means of inverse mapping 

M and P k 1 anH? <V ° r P rocedure in 3 ' D takes P la ce by the appropriate reordering of the indices' 

i, j, and k and (x, y, and z) coordinates of the physical space. 

1 n K°? C f P I ° f Spl i!, tmg n0t ° nly breaks com P Iex 3 -° flow field calculation into effective successive 

puters and ffeat l^geTn 1?"^ ad . vanta ^ of u the vector Processing architectures of modem com- 
n^v ? n d lt , lg , 3 ' D data f ts b y employing the block/pencil data base approach for large and com- 
P - ata (e.g Lomax and Pulliam (ref. 16) and Deiwert and Rothmund (ref. 17)) where the topo- 
logical space is subdivided into blocks which interface with one another exactly. P 


4.3 Boundary Modifications 

This section explains the adaptation procedure near or on the boundaries of the subdomain R Here 

Zr?n7* ,W ° r es: ™ dif,Cadon of equation (23) for the bound Jes of T and 2)7e m ai„,o 

nance of the configuration of the initial grid lines near the boundaries. For definiteness consider the nrin- 
cip e marching directions described in the previous sections. For k = k s the torsion forces T t , are 
^ualto zero, because the (k s - l)st plane corresponds to an external plane of the rerion R E^- 

hneTwould l thC 

ary Redf77ti77»riH U 7 , , iS T' 1 "?' prescrve the S™ 1 line configuration on or near the bound- 
e y ' • .? tb !° ° f g^ d P oints alon g ^jjc-lines may affect the direction of the initial grid lines n v or 
h ne ; gbborhood of the boundaries of R. This may be undesirable for cases whS for ins Sice 
h Un fvZ, ° a lS interface between two patched blocks of an initial grid, or when it corresponds to a 

GeSarr: ,his prob,em by 
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For clarity of discussion below let r , g , and n denote unit vectors along a reference line, a grid 
line, and a line normal to the boundary, respectively, issued from a grid point on the ^ lines on or near 
the boundary. As defined before, the reference line is a weighted average of vectors g and n ; then we 
can write: 

r=p!g+P 2 n (25) 

where (3j and p2 are some user specified parameters on the interval [0,1], (J 2 = 1 - Pi- Control of the 
vector r affects the inclination of the angles which in mm control the distribution of grid points along 
the adapted line i^k- To further clarify this point, suppose that in a given problem the user has to 
preserve the initial grid direction of the rij^-lines near the boundary, say for several ^-lines, stepping 
along the j-index from j = js to j = js + n g . This n g , where n g > js is specified by the user and is a 
mean to define a desired neighborhood near the boundary. Then by a simple algorithm, such as 

Pl->P, + (l-Pl)( " g: ^ — ] (26) 

one finds that f^r j = j s (i.e., at the boundary) Pi is increased to unity shifting the emphasis to the^irec- 
tion of vector g , then for j = j s + n g> Pi retains its original value, and for j s < j < n g , pi and g a^e 
smoothly modified. This procedure allows preservations of the initial directions of the rij^-lines of g . 
At the same time the value of the torsion constant Xi is modified in a similar way making the adjustment 
of Pi more effective. 

This argument can be repeated for problems where the initial boundary-conf omjiing grid must be 
preserved by smooth modification of P 2 . This would shift the emphasis to vector n near the boundary 
and smoothly transmit this effect to some specified neighbor of the boundary so that the orthogonality of 
grid lines near the boundary would be preserved. The procedure of modifying A.i and r just described 
for the r^k-lines can be repeated for all j-lines near the boundary plane £k> k = ks when stepping 
along the k-index. This approach makes the entire algorithm for treating the boundary consistent with the 
self-adaptive concept. 

There are also two more concerns with the boundaries. One, when the subdomain R does not cor- 
respond to the entire flow field, is that after each adaptation on the £k-plane the grid distribution along the 
lines ^jk, j = jE may not match smoothly with those ^-lines outside the boundaries, specifically the 
^jk- lines with j > jE- This case can be cured by projecting the grid spacing of the fij^-lines* j = jE onto 
the ^k-lines, j = jE + 1 and so forth, or by using algorithms such as discussed previously to preserve 
grid line configuration when j -> jE and k -> k£. Another concern is in avoiding large differences in 
spacing at the boundaries, along the J^-lines; in particular between Asj.i and Asjj^ for i = is, and 
between Asj_i and Asyjc for i = i£- These differences can be modified by the magnitude of tension 
forces like the weighting functions w 1Jv k at both ends prior to adaptation according to the concept of 
equidistribution; that is, we let 


i,j,k A s i-l,j,k - C , at i - is 

(27) 

i,j,k ^®i+l ,j,k — C , at i — 1 e 

(28) 


where C is an average of the products Wj j,k Asi j,k along the current J^-line (or previously adapted 
line Sj-i jc)- Then these values of Wj^k at both ends are merged to alleviate the values of several weight 
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functions in the neighborhood of i = is and i = ip The related flltrnrithm h 

self-adaptive concept and will be effective upon rtques, by ^ 


4.4 


Example of a Solution-Adapted Grid 


plume How is considered. Complex flow sm, pm' ™’ 1° jT! field about a tw »nozzle afterbody 

supersonic How and underexpanded supersonic ias^inthe beCaUSe of interaction between external 
mixing shear layers, barrel shocks, and Mach disk etc Thk features consist of shocks, 

geometric complexities and multiple iet ^am inre™eno v c ' ‘i further com Pl>cated due to 3-D 
computed multinozzle plume Hows with aTS XS Tl Va f a ‘ aP ?- y and Feiere isen ("* 1 8) have 
formed with an upwind solver Can ideal Li We XL okes . s ? lv «i thetr numerical simulations are per- 

initia. grid used b P y these authors^oSta ^h^nl C pSw't^'^ '° ^ 
under such free-stream conditions as- Mach numher - a 71 t 1 pl " ' Th com P utatl on is made 

atm. Figure 1 depicts the 3-D inirial grid £ the XL f“, : 2 ? K and press “ re = 0.0028 
the nozzle's exit plane. The How is bfTvmmetric ahnm th? ' d dle plume re ff lon downstream of 

solution in tennsof the Mach number. The Mach 

entire region; this means the function'^ eqtmtion (1 DhT^ 7 ,he sol “ don of Ma ch numbers in the 
of dte Mach number a!ong the direction^ Sradien, 

f *? ad T a,i ° n by the k 

boundary at k s = 47 k-index from the omj 

>index eS °" " * * P “ ~ by 

A. = oXTdX- aTOlTTpStX^On XfLTT “"r 'o' 5 *'* Asm1n = 0001 • As MAX = 1 , 
the initial grid line direction that is for the E '7dm "7° °[ R v ^ have requested the preservation of 
lows the specifications of the algorithm (e g^clusterino^H 1 ^ ^ ° veraI1 ada P ted grid clearly fol- 
barrel shocks, etc.) depicting "" Sh0CkS ' mixi " g la >™' 


4.5 


Example of an Initial Grid Generation 


an effc,CT„Sg7“d7ne™l P sT“ « of ' «• Captive gnd method as 

are related to the specific geoCCof C ^furtoTs^ S’”™" 0 "” 1 ^ wei S hIi "S fu"utions which 
sented here describes this procedure for a simnlp h vu 30 0 ri° w properties. The example pre- 
geometries can be subdivided into patched bliks of sfmnlp ^ 0metry ‘ sho u uld noted tha t complicated 
may separately be applied on each block. P g metry; then the following simple procedure 
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The surface geometry of the bump is defined by 


z(x,y) = < 


[1 + cos(4;ty)][l + sin(47t(x - 0.5) - 0.57t] , 

for 0.5 < x < 1 and 0 < y < 0.25 


( 29 ) 


0 , 


for x < 0.5 or 1 < x or 0.25 < y 


Initially a rectangular grid (49 x 29 x 40) is generated with uniform spacing on each side of a rectangular 
box B of the physical space; this is shown in figure 5. The coordinates T], and £ are now considered 
to be parallel with the direction of x, y, and z-axes. Next, the bottom surface of B is deformed to 
describe the bump, defined by equation (29), by shortening the z-coordinate accordingly. Figure 6 

rh of the grid construction in only one-half of the space. The computational domain R is 

ndices is = 1, 1 e = 49, js = 1 Je = 29, ks = 1, and k£ = 40. 


from figure 6, the % and r\ coordinates have conformed to the body shape, but have not 
the physical boundaries; £ coordinates along the z-axis are not orthogonal to the surface of 
make this grid more suitable for a 3-D flow field type of calculations, we need to properly 
; grid points in the physical space. In this instance, for a flow past a bump, we may cluster 
near the yz-symmetry plane along the T) -coordinates, and near the surface of the bump 
ordinates and adjust the inclinations of the £-lines so that they become orthogonal to the 
Dump. All this can be done by successive application of the described scheme as follows: 
the adaptation along the ^-lines on the ^-planes. Cluster the grid points using a 
F geometric weighting functions for regions of large curvature and exponentially stretched 
the boundaries. In equation (1 1) fy^ is replaced by a measure of the curvature and Mo 
in exponential function given by 


^i,j,k z i-l,j,k 2z i_i ,j,k z i+l,j,k » — 


ASt 


As 


ii 


(k/m) 


(30) 


-is - 1, Asi and Asp, the mesh spacings at both ends, are also specified. 

line inclinations are controlled by using large relative values for the torsion constants Xi 
; the grid quasi-orthogonal near the surface of the bump. 

id application of the scheme proceeds as previously mentioned schemes except for the 
g the riij-lines on the £-plane. The third application proceeds along the Cj^-lines on the 
fijjc in equation (30) set to zero because there is no curvature change along this direction, 
the final enhanced initial grid, which includes the features of clustering and orthogonality 
id boundaries of the bump. 


5. SUMMARY 


ation-adaptive grid based on variational principles and spring analogies was described. A 
s deal splitting procedure for multidirectional adaptation is used. Grid skewness and 

controlled by one-sided torsion springs. User-specified maximum and minimum grid 
spacings determine other important constants of the method and control the clustering and stretching of the 
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grid Tw ° computed examples were used to confirm the feasibility of the method and to demonstrate the 
adapuve gnd as robust and solution effective as well as a simplified initial grid generation scTeme 
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Figure 2 - 3-D initial grid-plume region. 
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Figure 4.- 3-D solution-adapted grid, based on Mach numbers. 
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Figure 7. Enhanced 3-D initial grid, based on curvature and stretch functions. 
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